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Abstract. This article solves an optimal control problem arising in attitude 
control of a spacecraft under state and control constraints. We first derive the 
discrete-time attitude dynamics by employing discrete mechanics. The orien¬ 
tation transfer, with initial and final values of the orientation and momentum 
and the time duration being specified, is posed as an energy optimal control 
problem in discrete-time subject to momentum and control constraints. Using 
variational analysis directly on the Lie group SO(3), we derive first order neces¬ 
sary conditions for optimality that leads to a constrained two point boundary 
value problem. This two point boundary value problem is solved via a novel 
multiple shooting technique that employs a root finding Newton algorithm. 
Robustness of the multiple shooting technique is demonstrated through a few 
representative numerical experiments. 


1. Introduction 

Typical space applications require reorienting a spacecraft or a satellite from 
a given initial configuration to a given final configuration over a given duration 
of time. Examples of such manoeuvres include positioning star sensors (attitude 
estimation sensors) towards deep space, pointing a camera in the desired direc¬ 
tion for imaging purposes, positioning solar panels for effective tracking of the sun 
for optimal energy harvesting, etc [34], Such orientation maneuvers are popularly 
known as attitude maneuvers, and the class of attitude maneuvers that optimize a 
certain performance objective are termed optimal attitude maneuvers. In most ap¬ 
plications, attitude control problems are subjected to additional constraints due to 
actuator saturation and/or forbidden regions of the state space. Control constraints 
are omnipresent — for instance, a momentum delivering device fitted on-board a 
spacecraft, e.g., a thruster, or a reaction wheel, has limitations in terms of max¬ 
imum torque produced and maximum momentum delivered. In case a spacecraft 
is fitted with a flexible structure like solar panels, it must inevitably have a re¬ 
striction of its maximum momentum for the protection of such equipments; such 
restrictions translate to state constraints in an optimal control problem. The class 
of constrained control problems involving orientation manoeuvres provides a rich 
supply of interesting and complex control problems related to space applications. 
Few computationally tractable solutions to such problems are known today, and 
in this article we propose a new and promising technique to solve a class of such 
constrained optimal control problems. 

In general, computationally tractable solutions to optimal control problems sub¬ 
ject to nonlinear dynamics are difficult to arrive at [35]. Constrained optimal at¬ 
titude control problems are vastly more challenging than standard optimal control 
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problems on Euclidean spaces because (a) the attitude kinematic equations evolve 
on the manifold SO(3), and classical techniques developed specifically for dynamics 
on Euclidean spaces, therefore, no longer directly apply, and (b) the presence of 
state and control constraints in the optimal control problem typically leads to two 
point boundary value problems (which represents first order necessary optimality 
conditions,) subject to inequality constraints, and these problems cannot be solved 
using conventional indirect techniques such as shooting methods [35] . In this article 
we provide an algorithmic solution to one such class of constrained optimal attitude 
control problems subject to state and control constraints, where the attitude kine¬ 
matics in discrete time are treated directly on the Lie group SO(3) and are derived 
using discrete mechanics [26]. 

Constrained optimal attitude control problems treated in the literature typically 
fall into one of two categories: one that considers only control constraints and the 
other that considers state constraints in addition to control constraints. Attitude 
control problems with control constraints have been discussed under the framework 
of optimal control with various performance indices involving minimization of time, 
fuel, and energy. A considerable body of work on the time-optimal attitude control 
problems exist [34, 29, 36, 10]. Fuel or energy optimal maneuvers are also of great 
interest since on-board energy sources are limited and precious, and there is there¬ 
fore a natural necessity for using them in an optimal way [3] . Early works on energy 
optimal control [2, 17] and fuel optimal control [9, 18] tackled the continuous-time 
optimal attitude control problem, arriving at first order necessary conditions for 
optimality via the Pontryagin maximum principle. The resulting boundary value 
problems were typically solved using shooting or neighborhood extremal methods, 
both of which suffer from high sensitivity to initial data. In recent years, attitude 
control problems with joint state and control constraints have been attacked from 
various directions. Attitude manoeuvres with joint state and control constraints, 
considering a linearized dynamics model is addressed in [37, 13]. In these works 
the problem is recast as a nonlinear optimization problem, and solved using direct 
multiple shooting methods or sequential quadratic programming [6]. Constrained 
attitude control problems with attitude kinematics represented in quaternion form 
was explored in [31, 24] using pseudospectral methods, and in [38] using particle 
swarm optimization. Representation of the attitude kinematics in quaternion form 
has a serious disadvantage of non-uniqueness; consequently, boundary conditions 
defined for attitude kinematics do not have a unique representation in quaternions, 
which presents a problem during computations. Another commonly used attitude 
representation is Euler angles, which suffers from the defect of singularity [33]. In 
order to avoid issues such as singularity and non-uniqueness, geometric techniques 
have been developed to solve constrained attitude control problems [26, 5, 21]. One 
such approach to attitude control problems in continuous time with control con¬ 
straints is discussed in [32]. Projection operators were used there to find a search 
direction on the underlying Lie group, while employing Newton’s method to solve 
the continuous time two point boundary value problem obtained via the Pontrya¬ 
gin maximum principle on manifolds [1, p. 165]. A different geometric technique 
to attitude control problem with state and control constraints is addressed in [14] 
where the authors attempted to handle state inequality constraints using penalty 
functions. As is well known, penalty functions do not enforce the state inequality 
constraints, but add penalties to the cost functions if the constraints are violated. 
This constrained optimal control problem leads to a two point boundary value prob¬ 
lem that represents first order necessary optimality condition. The boundary value 
problem is then solved by employing indirect single shooting method. None of the 
these works discussed the issue of digital implementation, a key aspect of which is 
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the process of discretization of the boundary value problems in time: Euler’s step 
and its derivatives are insufficient since these do not account for nor respect the 
underlying manifold structure of SO(3). 

In this article we tackle the issue of discretization up front by deriving a dis¬ 
crete time model via discrete mechanics [21, 19], and then employing variational 
analysis to arrive at first order necessary conditions for optimality. Discrete time 
models obtained via discrete mechanics are more accurate than other standard dis¬ 
cretization schemes such as Euler’s step because they preserve certain invariance 
properties like kinetic energy, momentum, etc, of the system, and the computations 
can be done directly on the manifold SO(3), (because this discretization respects 
the manifold SO(3),) thereby eliminating the problems associated with parametric 
representations of attitude. The presence of state inequality constraints in our at¬ 
titude control problem makes it challenging because the resulting boundary value 
problems obtained using variational analysis are subject to inequality constraints on 
the states. Such constrained boundary value problems cannot in general be solved 
using classical multiple shooting methods precisely because of those inequality con¬ 
straints. A non-classical multiple shooting technique has recently been proposed 
in [11] for solving constrained boundary value problems arising in optimal control 
problems with state constraints. There the complementary slackness conditions 
arising due to inequality constraints on the states are represented in the form of 
equality constraints using the Fischer-Burmeister function [11, Equation (2.7)] — 
a technique that can be highly inefficient in terms of computation because all the 
inequality constraints are considered at each iteration irrespective of them being 
active or not. In contrast, we propose a multiple shooting algorithm that deals with 
state inequality constraints in such a way that the dimension of the boundary value 
problem remains the same even when the state inequality constraints are active, 
making it more efficient in terms of time and memory complexity. 

We reiterate that the thrust of our contribution is towards computational trac- 
tability of optimal attitude manoeuvres. The advantages of our technique are five 
fold: First, the discrete time model derived via discrete mechanics avoids the need 
for discretization at later stages, and eliminates issues associated with parametric 
representations of the attitude kinematic. Second, an indirect multiple shooting 
method, which provides more accurate solution to the problem unlike direct tech¬ 
niques, is employed to solve the optimal control problem. Integration steps of 
the discrete kinematic equations represented by rotation matrices are calculated 
at each discrete instant, which provide exact attitude trajectories for manoeuvres 
unlike other standard schemes [35]. Third, the multiple shooting method that we 
employ here is more robust to initial guesses as compared to indirect single shoot¬ 
ing methods that are comparable in terms of accuracy [35]. Fourth, the proposed 
algorithm, being a multiple shooting method, can be implemented on a parallel 
architecture for fast computation. Fifth, the discrete time model obtained using 
discrete mechanics results in a boundary value problem that can be reduced to 
difference equations with momentum and co-momentum dynamics only. The di¬ 
mension of this reduced difference equation model is half of the original one, which 
further contributes to savings in terms of time and memory. Moreover, scaling of 
the co-state variables such that the co-states are invariant to step length selected 
for a given manoeuvre, improves the radius of convergence of the multiple shooting 
algorithm. 

This article unfolds as follows: In §2 we employ discrete mechanics to derive a 
discrete time model of the attitude dynamics. In §3 the energy optimal control 
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problem is posed as a discrete time optimal control problem, and first order neces¬ 
sary conditions are obtained using variational analysis. Then we discuss scaling of 
the variables and reduction of the dynamics to momentum and co-momentum vari¬ 
ables in §4. §5 contains an introduction to multiple shooting methods, and provides 
a solution to the system of difference equations represented in momentum variables 
presented in §4. §6 provides the numerical experiments for large angle manoeuvres 
with momentum and control constraints. We conclude in §7 with a brief discussion 
of future work. The proofs of our results are presented in a consolidated fashion in 
the Appendices. 

2. Derivation of Discrete time model via discrete mechanics 

This section contains the modeling of the attitude dynamics of the spacecraft 
using discrete mechanics [26] . In most optimal control problems involving mechan¬ 
ical systems, some sort of discretization is performed in order to employ numerical 
techniques. In the case of discrete mechanics, the variational description is directly 
discretized, and discrete time equations are obtained. This approach is advanta¬ 
geous in comparison to usual discretizations of continuous time models because it 
preserves certain invariants of the system such as momentum and energy. 

For ease of understanding, a quick introduction to discrete mechanics is given 
here, followed by discrete time modeling of the attitude dynamics of the spacecraft. 

2.1. Introduction. Consider a mechanical system with the configuration space Q 
as a smooth manifold. Then the velocity vectors lie on the tangent bundle TQ of the 
manifold Q and the Lagrangian for the system can be defined as L : TQ —>■ K. [25]. 
In discrete mechanics, the velocity phase space TQ is replaced hy Q x Q which is 
locally isomorphic to TQ. Let us consider an integral curve q{t) in the configuration 
space such that g(0) = qo and q{h) = qi, where h represents the integration step. 
Then, the discrete Lagrangian Ld '■ Q x Q ^ which is an approximation of the 
action integral along the integral curve segment between go and gi, can be dehned 
as [27] 

(2.1) Ld{qo,qi)~ ( L {q{t),q{t)) dt. 

Jo 

Pick h > 0, (this quantity plays the role of step length) and consider a grid of the 
time domain T = Nh as {tk = kh\k = 0,1,..., A^} and the corresponding discrete 
path space Vd{Q) ■= {qd ■ {tk}k=o Q}- discrete trajectory qd G Vd{Q) is 
such that qd(tk) = gfe. Now, defining the discrete action sum as 

N-l 

®d{qd) ■■= ^ Ld (gfc,gfe+i). 

k=0 

Assuming gd(0) = go and qd{tN) = qN fixed, define the variations S{qd) as S{qd(tk)) = 
S{qk) G Tq^Q which vanishes at the end points (i.e., S{qo) = S{qN) = 0). A discrete 
path qd is a stationary point of the discrete action sum ©^ if 'Dq^<3d{qd)d{qd) = 0 
for all 5{qd) [16, p. 21]. This is equivalent to saying that the points {g^} of the path 
qd satisfies the discrete Euler-Lagrange equations, i.e., 

(2.2) D 2 Ld{qk-i,qk) +DiLd(qk,qk+i) = ^ for all fc = 1, 2,..., A^ - 1, 

where Di is the derivative of the function with respect to the ith argument. 

Notice that (2.2) involves qk-i,qk and g^+i at A:th instant of time, which means 
that the difference equations obtained are of second order. To arrive at the discrete 
time model with first order difference equations, one needs to use the discrete time 
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analogue of the Hamiltonian formulation. To this end, recall that the continuous 
time Legendre transform is a map FL from the Lagrangian state space TQ to the 
Hamiltonian phase space T*Q. Similarly, the discrete time Legendre transforms 
¥'^Ld,F~Ld : Q X Q ^ T*Q [26] can be defined as 

F'^'Ld {qk,qk+i) {qk+i,Pk+i) = (gfc, £> 2^(1 (gfc, gfc+i)), 

F~Ld {qk,qk+i) {qk,Pk) = {qk, -DiLd {qk,qk+i )), 

which are maps from the discrete Lagrangian state space Q x Q to the discrete 
Hamiltonian phase space T*Q. The map F+L^ is the forward discrete Legendre 
transform which relates {qk,qk+i) to and the map F~Ld is the backward 

discrete Legendre transform which relates {qk,qk+i) to T*^Q. Let the discrete 
Lagrangian map Fl^ : Q x Q ^ Q x Q he defined as {qk,qk+i) = (9fe+i, 9fe+2); 
it defines the evolution of the dynamics on the discrete state space. Then the 
corresponding discrete Hamiltonian map : T*Q —>■ T*Q can be defined in the 
following equivalent ways [26]: 

Fl, := F^Ld o Fl, o (F±Ld)”' = F+L^ o {F-Ld )~"; 
this is clear from the commuting diagram shown in Figure 1. The discrete Hamil- 



Figure 1. Flow of the discrete Lagrangian and Hamiltonian map 


tonian map is defined in coordinates as follows: 
(2.3) 

{qk,Pk) Fl^ iqk,Pk) ■= iqk+i,Pk+i) where 


Pk = -DiLd (gfe,gfc+i), 
Pk+i = D 2 Ld {qk, qk+i) ■ 


2.2. Attitude dynamics in discrete time. We now apply the ideas introduced 
in §2.1 to obtain the discrete equations of the attitude dynamics of a spacecraft. 
First, the Lagrangian in continuous time is described, and then an approxima¬ 
tion of the continuous time Lagrangian is taken to define the discrete Lagrangian 
(2.1). Thereafter, the discrete time attitude dynamics is obtained using discrete 
Hamiltonian formulation (2.3). 

Consider a rigid body with a point, typically chosen to be center of mass, fixed 
on it. In order to define the orientation of a rigid body, two coordinate systems are 
considered with the origin at that fixed point. One frame fixed to the rigid body 
is known as the body frame, and the other is a frame fixed in space, known as the 
spatial frame. Let X be the position of the mass element in the body frame. Then 
the position of the mass element in reference frame x is related to the body frame 
coordinates X by the rotation matrix R{t) G SO(3) as x{t) = R{t)X. Let B be the 
region occupied by the body in its reference frame. Let p{X) be the density of the 
rigid body in the body coordinates at point X. Then the kinetic energy of the rigid 
body is [16, p. 243]: 
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which can be rewritten, in view of the left-invariance of the kinetic energy [7, p. 
275], as 


pm 


RX 


"d^X = ^J^p{X) 


R-^RX 


(RX. 


We know that the spatial angular velocity vector VL can be represented in terms of 
the body angular velocity vector w as 17 = R~^!jJ. Then 

(2.4) 17 X X = R-^uj X R-^x = R-^ (uj x x) = R-^x = R-^RX. 


Let 9 17 I—>■ 17 G so(3) be a vector space homeomorphism. Then, from (2.4) we 
arrive at the kinematic relation 17 = R~^R G so(3), the Lie algebra of SO(3). So, 
the kinetic energy can be represented in terms of the spatial frame angular velocity 
as 

K p{X) tr ((17X)(17X)^) (f’X = ^ tr (l7Jdl7^) , 

where 

Jd = i / p{X)XX^(fX. 

^ Jb 

The body moment of inertia matrix J ■= ^ p{X)X^XcPX is related to Jd by 
the following equation [23]: 

J = tr {Jd) /sxs — Jd- 

If the dissipative and potential forces are absent, then the Lagrangian L : TSO(3) —>■ 
IR, for the system is given by [16, p. 245] 

L(i?, 17) := R = i tr (l7 Jdl7^) . 

By the kinematic relation R = i?17 given above, we know that 17 = R^R. So, the 
Lagrangian can be written as 

(2.5) L{R, -R) = ^ tr (^R^RJdR^R) ■ 

We now proceed to discretize the Lagrangian (2.5). Considering discrete time in¬ 
stants tk = kh for k = 0 , 1 ,..., such that R{tk) = Rk and the approximation 
R{tk) ~ {Rk+i-Rk) ^ g \tk^tk+i] [23], the discrete Lagrangian Ld : SO(3) x 
SO(3) —IR is defined as: 

Ld{Rk, Rk+i) - hL (^Rk, 

h ( R^ {Rk+l ~ Rk) j {Rk+l ~ Rk)^Rk \ 

= 2 “■ (- h -*-5-j 

= ^ tr (^{rJ Rk+l - hy.3)Jd{R+ Rk+l - hxs)) 

( 2 . 6 ) =Ur{{l3^3-Fk)Jd), 

h 

where Fk = R^Rk+i- Note that under the discretization technique employed, the 
discrete Lagrangian, like its continuous counterpart, is invariant under the action of 
the SO(3) group. This property will be useful later when momentum equations will 
be derived and the rotation sequence constructed based on the momentum history. 

Our objective is to come up with first order difference equations describing the 
attitude dynamics of the spacecraft. By the left trivialization of the cotangent 
bundle of a Lie group, T*SO(3) can be represented as SO(3) x so(3)*, where so(3)* 
denotes the dual of the Lie algebra so(3) [1, p. 254]. We now proceed to find the 
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discrete time Hamiltonian map (2.3) : SO(3) x so(3)* ^ SO(3) x so(3)* such 

that 

f Hfc = —DiLd {Rk, Rk+i), 

Fl^ = (^i?fc+i,nfe+i) , where < n^+i = £>2^^ {Rk,Rk+i), 


In order to find H^ and H^+i, the variations in Rk are defined in terms oirjk € so (3), 
and the expressions DiLd {Rk, Rk+i) and D 2 Ld {Rk, Rk+i) are evaluated using the 
duality product on so(3) [25, p. 290]. 

For a given e S IR and r]k S the variation in Rk can be defined as 


(2.7) 


5Rk '■= -j- 
ae 


Rk = Rkfjk, where rjk e so(3); 


e=0 


then the duality product of H*, e so(3)* and rjk S so(3) is defined as 


tr (fi-kfij^ =: (fik,Vk^ 


de 


1 


Ld{Rl,Rk+i) = - tr (((5i?fc i?fc+i) Jd) 


e=0 

= {dlRjRk+iJd) = ^tT (FkJdvJ) ■ 


Hence, 


tr 


V 


iflfc - ^FkJd ]vk =0 for all rjk G so(3), 


Ck 




which means that Ck is a symmetric matrix. Its skew-symmetric part is, therefore 
zero, and this leads to 


( 2 . 8 ) 


hUk = FkJd — JdFk ■ 


Similarly, the duality product of H^+i £ so(3)* and 77^+1 G so(3) gives. 


tr =; (^Hfc+i, 77 ^+ 1 ^ = 


de 


1 


Ld{Rk,Rl+i) = rfo (- {Rj6Rk+i) Jd) 


e=0 


and 


tr 


= (i?^i?fc+i 77 fe+iJd) = i tr (JdFfc? 7 fe+i) , 


Hfe+i - ^JdFk ] 77^+1 = 0 for all 77fc+i G so (3) 


V 


Dk 




which means that Dk is a symmetric matrix with its skew-symmetric part equal to 
zero. Therefore, 


(2.9) Hfc+i = = FjlikFk = FflTfe, 

leading to the following update equation for the momentum: 

( 2 . 10 ) Uk+i = Fjnk. 

In the presence of control, (2.10) modifies to 

( 2 . 11 ) Uk+i = Ff[llk + huk. 
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where Uk is the control input at fcth instant of time. The rigid body equations in 
discrete time are finally obtained as: 


Rigid Body 
Dynamics 


{ R/c+l— 

Bfc+i— Il/g “t” huf^^ 
hUk — FkJd — JdFjF. 


3. Optimal control of discrete time attitude dynamics 


We state the optimal control problem arising in executing energy optimal atti¬ 
tude manoeuvres of a spacecraft. The spacecraft is assumed to have three actuators, 
aligned along the three principal moment of inertia axes. Each actuator has its indi¬ 
vidual saturation limits. The objective is to find the energy optimal control profile 
for orienting the spacecraft from an initial configuration to a desired configuration 
in a given duration of time, while obeying pre-specified momentum bounds. First 
we pose this requirement as an optimal control problem, and then derive the first 
order necessary optimality conditions using variational analysis. Later §5.2 the 
boundary value problem obtained as the first order necessary conditions will be 
solved using a novel multiple shooting method. 


3.1. Problem description. Our objective is to find the energy optimal control 
law to manoeuvre a spacecraft from the initial configuration (i?i,ni) to the final 
configuration {Rf,Ilf) in N discrete time steps satisfying the following constraints: 


(1) < F k = 0,1,..., N — 1, and i = l,2,3, 

(2) n^|<d* k = l,2,...,N-l, and z = l,2,3. 

This problem can be posed as an optimal control problem in discrete time as follows: 


(3.1) 


1 

minirnized V 


subject to 

(3.2) 

f Rk+1= RkFk 

system of equations < nfc+i= F^ 11^ -|- huk with 

[hllfc = FkJd — JdF^ 


(3.3) 

boundary conditions (i?o, IIo) = (i?i, 11^), (R^r, IIat) = (i?/, 11/), and 

(3.4) 

i (ul)'^ < (F)'^ for all fc = 0,1,..., fV — 1, and *=1,2,3, 
constraints < ; ,2 1,2 

[(n^) ^ i^’’) for all fc = 1, 2,..., iV — 1, and *=1,2,3. 


Note that the optimal control problem (3.1) has both control and state inequality 
constraints. While the individual control inputs are constrained in magnitude, the 
performance measure reflects a 2-norm on the control action at each stage. 


3.2. Necessary optimality conditions. We represent the variations of Fk in 
terms of variations in 11^ and then the first order necessary conditions are derived. 

• Representation of the variations: 

Using (2.7) the variations for the matrix RjRk+i is defined as 

^ i^Rk Rk+1^ — ^Rk Rk+1 F Rf^ — VkRk T RkVk-\-l- 


(3.5) 
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Using the property x = F^xF, (3.5) simplifies to 
(5 (R^Rk+i) = Fk 77fc + 'nk+i)'' , 

where (■ )^ : so(3). Similarly, for a given S we define the 

variation in Fk S SO (3) as 

^Fk = Fkik- 

The implicit equation hX\k = FkJd — JdFk in (3.2) gives the relation be¬ 
tween momentum and change in orientation at fcth time instant. So, the 
relation between the variations in momentum (511^ and 5Fk can be obtained 
from the implicit equation as [21] 

h5Wk = SFkJd — JdSFk = FkikFkJd + JdF^ Fk^k, 

which can be further simplified using xA -|- Al^x = ({tr[y4]/3x3 — A} x)^ to 

Mllfe = ((tr [FkJd) hx 3 — FkJd) Fkik)'^ ■ 


Lemma 3.1. The matrix (tr (FkJd) Isxa ~ FkJd) is invertible if 

, /2d3 + d2-di 

V 2 j ^ V 2(d3 + d2) ’ 

where Fk = e^'‘, ^k G R^, and Jd = diag(di,d 2 ,d 3 ) such thatO < di < 

c?2 < da. 

We present a proof of Lemma 3.1 in Appendix A. Armed with Lemma 
3.1, we represent the vector fk in terms of the variations in momentum i.e. 
SUk [21] as 

fk = BkSUk, 

where Bk = hF^ (tr (FkJd) hxa - FkJd)~^ ■ 

• Necessary optimality conditions: 

Let Xk € R^ and \k G R^ be the Lagrange multipliers corresponding to the 
equality constraints i?fc+i — RkFk = 0 and Ilfc+i — F^Ak — huk = 0 respec¬ 
tively. Similarly let 0 < £ R and 0 < £ R be the Lagrange multipliers 

corresponding to the inequality constraints < (c®)^ and < 

(d*)^. Let us justify why the Lagrange multiplier Xk £ R^ is chosen cor¬ 
responding to the rotational kinematics Rk+i — RkFk = 0. Rotational 
kinematics can be rewritten as RilRk+i — Fk = 0, where Fk can be identi¬ 
fied by its skew symmetric part, which in turn can identified by a vector in 
R3 [22]. 

Claim 3.2. Consider the equality RilRk+i = Fk- If we assume that the 
step length h is small enough such that the relative orientation RilRk+i 
between two adjacent time instances tk and tk+i is less than i.e., 

Il^fcll) IlCfell < ^ where RJRk+i = ef^ = e^'= and fk,fk £ R^, 

then the aforementioned equality is satisfied if and only if the skew sym¬ 
metric parts of both sides are identical. 
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The augmented performance index can be defined as 


minimize Za 


:= — (ukiUk) + (Afc, —Hfc+i + Ilfc + huk) 

k^O 

+ (xk, i {Fk - {RjRk+i - Rj+^Rk)"'^ 

1 1 
+ 2 i^k, Uk&Uk-cQc) + ^ - (/3i, 11; 0 IIj - d 0 d), 


where 

c = (c^ d={d} (P , xQy=(x^y^ x'^y'^ x^y^)^. 

Using (3.6), (3.7) and (3.8), the infinitesimal variation of the augmented 
performance index is defined as 


AT-l 

<5aa ■.= '£{{Suk,Uk) + (^^,6{Fk-F^y -d{RjRk+i-Ri+,Rkf) 

k^O 

+ (Afc, —+ SFjI Ilfc + F^SHk + hSuk) + {Suk^ ak (D Uk)} 

N-l 

fY.{SUi,PiQUi), 

1^1 

where 

(3.9) al,l3l>0 for all i = 1,2,3, 
together with the complementary slackness conditions 

(3.10) /3? ((ni)2 - (d*)2) = 0 and ^ ((4)^ - (c^) = 0. 

Employing the property that xA + A^x = {(tr[A]/ 3 x 3 — A) a;}^ and that 
the variation on the boundary is zero, i.e., rjo = = 0 , dllo = dllAr = 0 , we 

rearrange the terms of the expression 5Za and get the following expression 
after standard algebraic manipulations: 

A-l 

dZa = F, {^'k^k,Uk + hXk + akQ Uk) 
fc=0 
A-1 

+ X! {Fk)hx^ - Fk)xk - (tr (Ffc_i) /3 x3 - Fk-i) Xk-i^ 

k^l 

+ ^ -Afc_i + (Fu - BlF^k) Xk) 

k^l 

N-l , . 

(3.11) + ^ (dnfc,,afe 0 nfe + - 6 ^ {t'^{Fk)h^z-Fk)xk) ■ 

k^l ^ ! 

By first order necessary condition of optimality, SZa = 0 along all possible 
variations (5zifc, 77 ^, 511^, and we obtain the co-state equations from (3.11) 
as: 

(3.12) 

Co-state j i^k-l) ^Sx3 t^k—l) Xk—l ^k (tr (-f/c) ^3x3 ^k) Xk — 0: 

Dynamics 0 - Xk-i + - BJfFF) Xk + \BI (tr (F^) / 3 X 3 - Fk)xk = 0, 
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and the optimality condition for the control as 

(3.13) Mfc +/lAfe + afe © Mfe = 0. 

From (3.10) we know that if u\ lies in the interior of the feasible region 
C* := {ii® G 111(11®)^ < (c®)^}, then a\ = 0. So, by (3.13) we have hX). = 
—Uk- On the other hand if u\ = c* then by (3.13) and (3.9) we have hX], = 
—c*(l + a].) < —c\ Similarly for = —d, we have hX\ = c*(l + a].) > d. 
Hence, the optimal control can be written in a compact form as 

(3.14) =-min{c\ |/iA^|}sgn(A^). 

We represent the state and co-state dynamics in terms on momentum and co¬ 
state corresponding to momentum variables. This technique will reduce the model 
of the system and hence the algorithm will perform better in terms of memory and 
time requirement. 


4. Scaling and model reduction 

The orientation boundary constraints (Rq^Rn) = can be represented 

in terms of Fk, which in turn can be computed for a given H^ using the implicit 
form 

(4.1) mk = FkJd - JdF^ ■ 

We represent orientation constraints first in terms of Fk and then we discuss a 
technique for computing Fk for a given value of Hfc. 

4.1. Representing orientation constraints in terms of momentum. We rep¬ 
resent the boundary constraints on orientation {Rq, Rn) = {Ri, Rf) (3.3) in terms 
of Fk- Then using (3.2) we reconstruct the orientation Rk variables, we see that 
Rf can be represented in terms of for A: = 1, 2,... , — 1, as 

Rf = Rn = RiFoFiF2 ... Fn-i. 

The boundary constraints of orientation can be rewritten as: 

(4.2) RjRf = F 0 F 1 F 2 ...FN- 1 . 

To represent the constraints (4.2) in vector form, we have nine equations. From 
(3.12) we know that the number of free variables corresponding to orientation kine¬ 
matics are actually three: xo G So, the boundary conditions (4.2) have to 
be represented by three independent constraints. Let us define the maps SO(3) 9 
M i-> logm(M) G so(3) and so(3) 9 a: >->■ {x^ G R^. Then the boundary condition 

(4.2) is satisfied if 

(4.3) Orientation constraints: Comt := (logm (i?Ji?iFoFiF 2 ... = 0. 

We need to compute the gradient of the orientation constraints (4.3) w.r.t. mo¬ 
mentum Hfc. This can be done as follows: From (4.3) we know that 

= Jij R^FqFi ...Fk... Fat-i, 

leading to 

e'^ornt = RjRiFoFi.. .X>n'^(^fc) ■ • 

After algebraic manipulations we get 

^ni^ornt = ^IV-l^JV-2 ' • ' (^fc)Ffc+i . . . F^-l, 

and using the property A^x = A^xA we conclude that 

RniCornt = F^_iF^_ 2 . . . (^FJVui^Fk'^ . 
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Therefore, the gradient of the orientation constraints (4.3) w.r.t. momentum Ilfc is 
(4.4) ^Hfe^ornt — (^n^^ornt ^H^^ornt ^n|^ornt-) 

We notice that the orientation constraints (4.3) can be computed only when we 
construct the matrix for a given value of the momentum vector Ilfc. Similarly, 
the gradient of the orientation constraints (4.4) can be computed only when we find 
the derivative of Fk w.r.t. for fc = 1, 2,..., iV — 1, and i = 1, 2,3. 


4.2. Determining Fk and V-ni Fk in terms of momentum Ilfc. In order to 
calculate the orientation constraints (4.3), Fk is obtained in terms of Ilfc. We know 
that Fk can be obtained from Ilfc by solving the implicit equation (4.1). To solve 
this implicit form, we choose quaternions to parameterize the matrix Fk. Let 

Ao + 9i - 92 - 2(7192 - 2(7 o93 29193 + 29092 \ 

Fk := 29192 + 29093 9o - 9i + 92 - 93 29293 - 29091 

\ 29193 - 29092 29293 + 29091 9o - 9i - 92 + 93 / 

and the inertia matrix J and Jd are defined as: 


/F 

0 



^_F + F + F 

0 


'll 

0 

F 

0 , 

0 

1—1 

1 

1—1 

+ 
1—1 

0 

\ 0 

0 

r 

2 \ 

V 0 

0 

F + F - F/ 


Let Ilfc := (n^ be the momentum of the body at the fcth instant, and 

9 = ( 90 ) 9i: 92 , 93 )^; then (4.1) can be represented in the form of nonlinear algebraic 
equations as follows: 


(4.5) 


9(9(nfe),nfc) 


/29293 (I" - I") + 290911" - ^ni\ 

29193 (F-F) + 29o92F-/in2 
29192 (F-F) + 29o93P-/in3 

\ 9o + 9? + 9i + 9i - 1 / 


= 0 . 


For a fixed value of momentum Ilfc, the system of equation (4.5) has quaternions 
as unknown parameters which can be numerically found using Newton’s method 
at each instant of time. Since Fk is represented in terms of quaternions, we first 
need to find the variations of the quaternions in terms of the momentum vector Ilfc; 
these can be obtained by taking derivative of (4.5) w.r.t. Ilfc as 


( 9 (nfe), Ilfc) T>nfc 9 (nfc) + ( 9 (nfc), Ilfc) = 0, 


resulting in 


/ dqo 


(4.6) 


T^qg (9(nfc),nfc) 


dqi 

3(?2 

dqs 

anF 


dqo 

dqi 

dq-i 

dqo 

dm 


dqo 

dqi 

dqo 

dqo 



fh 

0 



0 

h 

0 


0 

0 

h 


^0 

0 

0 / 


'Dn^qijlk) 


where 


'Dqg (9(nfc),nfc) 


/291F 

290F 

293 (F - F) 

292 (F - F)\ 

292F 

293 (F-F) 

290F 

291 (F-F) 

293P 

292 (F - F) 

to 

1 

1—1 

290F 

\ 290 

291 

292 

293 / 
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The matrix 'Dui^q can be obtained by solving the linear system (4.6). The derivative 
of the matrix Fk w.r.t. 11^ is obtained using the chain rule as 


(4.7) 


Vy 


Fk 


n—0 


dqn 

dUl 


We now discuss the reduction of the difference equation model (3.2), (3.12) to the 
momentum and comomentum dynamics respectively and later the reduced model 
is scaled by change of variables so as to make the model invariant under change 
in the step length h. Invariance of the difference equation model means that for 
a particular manoeuvre the optimal trajectory and corresponding Lagrange multi¬ 
pliers remain identical for different step lengths h. This matter is quite essential 
because it largely affects the region of convergence and order of convergence of the 
algorithm [12]. 


4.3. Scaling and model reduction. First we discuss about model reduction and 
later scaling of the reduced model by appropriate change of variable. Let us define 
a new variable Cfc := 2h“^ (tr {Fk) hxs — Fk) Xk then (3.12) can be written as 

(4.8) /3fc 0 Hfe - Afe_i + (Fu - hMkl^k) Afc + ^Cfc = 0 

(4.9) Cfc-i — FkCk = 0 

where Nk = From (4.9) we conclude Cfe = Qk Co such that Qk = F1F2 ■ ■ - Fk- 
So (4.8) and (3.2) can be further reduced to the system of difference equations 

(4.10) 

Ilfc+i Fj^ huk 

Pk+i © Ilfc+i — Afc -f ^Ffc+1 — hAfk+iFf^^^Ilk+i^ Afc+i -f A4-I-1 ^Ty^Co 

Assume for the moment that for a manoeuvre the fth element of the control vector 
Uk saturates (i.e., |u^| = F) at the fcth instant of time. Clearly, from (3.14) 
-^fc > length h is small, A^ will take very large values, which makes 

the difference equations (4.10) stiff. To avoid this situation, we define a variable 
7fc := /lAfe, and modify the difference equation (4.10) to 

(4.11) 

Ilfc^i F). Ilfc huk 

hPk-\-i © Iffc+i 7 fc “t“ ^Ffc_|_i /iAffc+iF^_j_^Ilfc_|_i^ Xk-\-i T -^k+iQk-i^iCo 

where 

(4.12) =-min{c\ |7(,|}sgn(7^), A4 = {ti {JdF,[) - JdF^) ^ Fk. 

Note that Fk and its gradient with respect to momentum can be obtained as dis¬ 
cussed in §4.2. In the following section we employ multiple shooting method to 
solve the system of difference equations (4.11) with boundary conditions (3.3) and 
constraints (3.4). 




5. Multiple Shooting Methods 

Shooting methods were mainly developed for solving ordinary differential equa¬ 
tions or difference equations with given boundary conditions. An initial guess is 
taken for the unknown initial values of the differential or difference equation vari¬ 
ables. Then the variables are computed at the terminal time and compared with a 
known value of the variables at the boundaries. Then the initial guess is improved 
at each iteration to match with the known boundary values. The multiple shooting 
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methods, the time domain is divided into sub-intervals (time domain decomposi¬ 
tion), and the boundary value problems are solved for each sub-interval with the 
condition that the boundary values at the common points of the adjacent intervals 
are the same. Multiple shooting methods are generalizations of the single shooting 
method in the sense that multiple shooting with a single time interval is equivalent 
to the shooting method [15]. 

Multiple shooting methods have many advantages over single shooting meth¬ 
ods. The former are more stable, and hence can be applied to stiff problems, and 
the time domain decomposition allows one to introduce the initial guess to the 
problem with prior knowledge. Multiple shooting methods allow one to compute 
the solution of the differential equation at individual intervals, which can be very 
efficient in computation using parallel architecture [15]. 


5.1. A quick introduction to multiple shooting methods. Multiple shooting 
methods constitute generalizations of the single shooting method, in which two 
point boundary value problems are solved at each iteration for the subintervals of 
time domain simultaneously. Let 

(5.1) x = f{x,t) cc G R" with boundary conditions 

(5.2) R : R" X R" ^ R" such that B{x{a),x{b)) = 0. 

Let the time domain be decomposed into N sub-interval as 


to = To < Ti . . . Tat-i < Tn = tf, 


and let us consider {N + 1) variables s°, ..., s", known as the multiple shooting 

variables. These multiple shooting variables are the guessed initial values of the 
dependent variable x defined in (5.1) at the specified time instants. We now define 
initial value problems for each sub-interval as 

(5.3) = f{x^,t) such that : [Tfe,Tfe+i] R" with 

(5.4) x'^{Tk) = s^ forfc = 0,l,...,(A-l). 

Note that the solution {x^( •, s^)jA: = 0,1,..., {N — 1)} to the initial value prob¬ 
lems (5.3) can be a solution to the boundary value problem (5.1) only if the solution 
x^ of the interval [Tfc,Tfe+i] matches with the initial condition for the next inter¬ 
val i.e. x^(Tfe_|_i, s^) = This condition is known as the matching condition. 

These matching conditions for each interval can be combined together and can be 
represented in stacked form as 


(5.5) 


F{s) 


/ s^-x°(ri,s°) \ 

- X^{t2,S^) 


gN _ xN 1 ^ 

\ B{x{a),x{b)) j 


= 0 , 


with x^{ •, s^) the solution of initial value problem (5.3), and s := (s°, s^,..., s^). 
The algebraic equations (5.5) can be solved using Newton type algorithms for multi- 
variable functions is described in Algorithm 1. 
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Algorithm 1 Multiple shooting method 
1: procedure Root of the function F 

2: Choose appropriate initial guess s and the error tolerance bound d > 0. 

3: Stopping Criterion: 

4: if ||i^(s)||^ < (5 then Stop. 

5: end if 

6: Calculate As by solving the linear system 'DsF{s)As = —F{s). 

7: Update the value of s = s + As. 

8 : goto Stopping Criterion. 

9: end procedure 


5.2. Multiple shooting method for attitude dynamics. The system dynam¬ 
ics (4.11) with the boundary conditions (4.3) and (ftojllAr) = (11^,11/) have to be 
solved with additional inequality constraints (ft^)^ < (d*)^ along with the compli¬ 
mentary slackness conditions (3.10) defined by 

(5.6) 

( llfe+i - F/lUk - huk = 0, 

(5.7) 

boundary conditions |nAr — ft/ = 0 , Hq — Hi = 0 , Comt = 0 , 

(5.8) 

slackness conditions |/3^ ((H^)^ - (d*)^) = 0 , /3^ > 0 , (H^)^ < (d*)^, 

where 

III = - min {c\ \jI I } sgn( 7 ^), A4 = (tr {JdF ^) /gxs - JdF ,[) ^ Fk, 

Cornt := (logm [rJR iFoFiF 2 ... Fn-i))'^ ■ 

First we consider the case in which state constraints are not active (i.e., = 0 

for all k). Then (5.8) is trivially satisfied. If state constraints are not active, then 
the necessary conditions are defined by (5.6) with /3fe = 0 for all k, and (5.7). We 
now represent the matching and boundary condition for the system of difference 
equations (5.6) and (5.7) in terms of life, 7 ^ and Co and solve the system of nonlinear 
algebraic equations comprises of matching conditions. 

5.2.1. Matching conditions for multiple shooting methods. The system of difference 
equations (5.6) along with the boundary conditions (5.7) will result in the following 
set of matching conditions, which can be solved using Newton’s root finding algo¬ 
rithm with quadratic convergence rate [28]. 


Solve: 



(5.9) 

7W(A):=(Si Si ... Efc Sfe 

Eiv ■^N Cj]jt,j„ — 0 

where 

X-.= {U^ 7oT ... rT ... 

U/c+i .— (R/c-t-i fffc dzifc) , 

n^v 7iv (d) I 


Cmtm := I , Cornt := (logm {rJR,FoFiF 2 ... Fn-i))"' . 
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Figure 2. Newton’s root finding algorithm 


Let VxAi be the Jacobian matrix of the matching conditions Ai (5.9). Assuming 
Xn to be the solution of the system A4 at the nth iteration of the Newton’s root 
finding algorithm, the (n + l)th iteration is given by 

Xn+l = Xn + AA„, 

where AA„ is a solution of the linear system X>xAJ(A„)AA„ = — A1(A„). 

It is important to note that Newton’s update AA„ is not necessarily unique; 
AXn is unique if and only if 'DxXi{Xn) is invertible. Invertibility of the matrix 
'DxM is proved for a special case (when only momentum dynamics are considered) 
in Appendix C. 

Notice that Uk defined in (4.12) is not differentiable, but it is Lipschitz continu¬ 
ous. So, we take its generalized gradient [8], defined by 

f-1 if7fc < 

= S -1 if -7fc > -c^ fori = 1,2,3. 

[ 0 otherwise 

The resulting nonlinear algebraic equations (5.9) can be solved using the non¬ 
smooth version of Newton’s method in [30] , and this is illustrated with the help of 
the flowchart shown in Figure 2. 

5.2.2. Multiple shooting method with state constraints. Here we discuss the general 
case in which momentum constraints are active. By the complementary slackness 
conditions (5.8), it is clear that if the momentum constraints are not active, then 
/?(, = 0 for all i = 1,2,3, and k = 1,2,...,A — 1. On the other hand, if the 
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slack variable > 0, then the momentum variable corresponding to the slack 
variable /3^ lies on the boundary (i.e., 11^ = ±d*) [6]. 

The original Newton’s root finding algorithm cannot handle momentum inequal¬ 
ity constraints; we propose a modified strategy for handling the momentum in¬ 
equality constraints. The algorithm has three phases: 

(1) Projection of the states to the feasible region: After each Newton’s 
iteration, the momentum variables {Ilfej^Q are projected onto the feasible 
domain. 

(2) Identification of active constraints: Active constraints can be identified 
by computing from (5.6). If > 0, then the inequality constraint 
corresponding to the state 11^ is active; else it is inactive. 

(3) Enforcing active constraints: When /3^ > 0, then the equation with 

in (5.6) is trivially satisfied. So, the matching condition E], will be replaced 
with the corresponding active constraint |n^| = d® for the next iteration. 

Let 


A:=(nJ ... nj ... nl 7^ Co^)^ e Reiv+g 

be a vector in the augmented space; the feasible region A is defined as 


A := {A G 1R®^+9 I |n*fc| <(f}. 


It is clear that A is the intersection of the two half spaces A+ := {A | 11^ < d* } 
and X~ := {A | 11^ > —d®}, hence a convex set. Let Va : 1R®^+® —A be the 
projection onto the feasible region A, and M represents the active constraints, 
which is a collection of equality state constraints, active inequality state constraints 
and co-state constraints corresponding to inactive inequality state constraints. Our 
modified multiple shooting method is described in Algorithm 2, and is illustrated 
with the help of the flowchart given in Figure 3. 


6. Numerical Experiments 


The following data have been considered for the simulations: Principal Mo¬ 
ment of inertia of the satellite (P,P,F)=(800,1200,1000) kgm^, Sampling time 
(T) = 0.1 s, Range of angles {9) can vary between [10°, 90°] about any axis, Max¬ 
imum torque or control bound (c) = (20, 20, 20) Nm, maximum momentum (d) = 
(70, 70, 70) Nms, time duration (tmax) can range from 5 s to 30 s. 


An open loop optimal control profile is obtained for a maneouvre of reorient- 


from initial momentum 11, = 


ing the satellite for 90° about axis 

(30, —10,10) Nms to desired momentum 11/ = (0,0,0) Nms in 19 s. The optimal 
control profile along with the momentum and co-state vectors is shown in Figure 4. 
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Algorithm 2 Modified multiple shooting method 


procedure Root of the function M such that A g A 

Choose appropriate initial guess X and the error tolerance bound <5 > 0. 

Projection: 


6 : 

7: 


Compute X i-A Vk{X) := (nj 7 ([ ... flj 

where fl^ = sgn(n^) min {d*, |II^ |}. 

Using (5.6), update co-states 


■ • ■ 


7fe-i = 


-ui if|M*fc|<c\ 
Ik-i elsewhere, 


where Uk = 




if ni = d® then 


/3fc = 


Compute the slack variable using (5.6) as 

^Z(nfc)E!/j i(nfe Co)) 




, where ^(11^) = 


Fk 


-MkFjnk]. 


end if 

Define active constraints 
7W(A):=(Si Si ... Efc S^ 
where 






'-'ornt/ 


“fc 



if/3^ >0, 
otherwise. 


10 : 


11 : 

12 : 

13: 

14: 

15: 

16: 


Evaluate the matching conditions M{X) and its gradient VxXi(X), 
whereA:=(nJ jT jjT ~t ^ _ jjT ~t 


if 


M(X) 


end if 


< (5 then Stop. 

OO 


Calculate AA by solving the linear system 'DxM(X)AX = —Ad(A). 
Update the value of A = A -I- A A. 

goto Projection, 
end procedure 





Figure 4. Rotation of 90° about axis with ini¬ 

tial momentum fli = (30,-10,10) Nms, desired momentum lly 
= (0, 0,0) Nms, time duration = 19 s. 
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Figure 3. Modified multiple shooting algorithm for problems 
with state constraints 


It is clear from the Figure 4 that the control saturates when the absolute values 
of the co-state vectors corresponding to the momentum are more than the control 
bounds. If the momentum saturates, then the control along that axis goes to zero. 


7. Conclusion and Future Direction 

A tractable solution to a class of optimal control of spacecraft attitude ma¬ 
neuvers under control and momentum constraints has been presented. We used 
discrete mechanics to discretize the continuous time model, which has many advan¬ 
tages over the conventional discretization schemes like Euler’s steps. In particular, 
the model so obtained preserves conserved quantities of the body like momentum 
and energy. It is efficient in terms of numerics because the model reduces to the 
momentum dynamics only. A new multiple shooting algorithm has been proposed 
to solve the system of difference equations obtained as first order necessary con¬ 
ditions using variational analysis. This algorithm can be used to solve optimal 
control problems with state inequality constraints. In future, convergence results 
for the modified multiple shooting methods will be explored. The proposed multi¬ 
ple shooting method uses Newton’s root finding algorithm at each step for finding 
the roots of the matching conditions. Newtons’s method has a quadratic rate of 
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convergence in a neighborhood of the true solution, and therefore it is worth ex¬ 
ploring the convergence rates of our algorithm in a neighborhood of the matching 
conditions. It is important to note that our algorithm can handle box constraints 
only; an extension to more general constraints is under development. 


Appendix A. Proof of Lemma 3.1 
Lemma 3.1. The matrix {tr (FkJd) hxa — F^Jd) is invertible if 

||^fc||\ , / 2d3 + d2 — di 


cos 


< 


2(^3 + 1 ^ 2 ) 


where Fk =6^*", fk € and Jd = diag{di,d 2 ,d 3 ) such thatO < di < d 2 < d^. 


Proof. If ti (FkJd) = 0 then the matrix {tr {F^Jd) hxs — FkJd) is invertible. If 
tr {FkJd) ^ 0 then {tr {FkJd) hxs — FkJd) is invertible if and only if the ma¬ 
trix ^/ 3 x 3 — trfpkJd ) ) invertible. Using Banach lemma [20, p. 193], the matrix 


(r., _ FkJd ] 
\^3x 3 tr(FkJd)) 


is invertible if 


FhJd 

tr{FkJd) 


< 1 . 


FkJd 
tr (FkJd) 


\tr(FkJd)\ ~ \tr{FkJd)\ 


Then 


FkJd 

tr(FkJd) 


< 1 


holds if 


||Jd|| < \tr{FkJd)\. Now we show that if cos < 


y tldf+d^)^ then condition ||Jd|| < |tr(FfcJd)| is true. Let us choose the quater¬ 
nions {qo,qi,q 2 ,q 3 ) G R'^ as a parametrization for the rotation matrix Fk where 
go = cos and (gi,g 2 ,g 3 )^ = CfcSin the vector represent the unit 

vector corresponding to vector fk- Then 


|tr {FkJd)\ = |di(go + <7? - 92 - Qs) + ^ 2(90 - 9? + 92 “ qI) + <^ 3(90 “ 9? - 92 + 93 ) 

> di(9o + ql-ql- ql) + Mqo -ql + ql- ql) + Mqo -ql-qlF ql) 

= di + d 2 + ds — 2 {(di -|- ^2)91 + (^2 + d 3 )q^ F {d^ F d.i)q^^ 

(A.l) > di -|- c?2 + 1^3 — 2/, 

where 


(A.2) / = maximize (di -I- d 2 )x -I- (d 2 -I- d 3 )y F [dz F d/fjz 

x,y,z 

subject to 

xFyFz=l-ql, 
a; > 0, 2 / > 0, z > 0. 


The optimization problem defined in (A.2) is a linear programming problem and 
hence the optimum value will be attained on the vertices of the feasible region. 
Without loss of generality, assume that 0 < di < d 2 < d^. Then / = (d 3 -|-d 2 )(l—9o). 

Given that cos 9o = cos one can conclude that 

^0 < ^tldf-hdF^ • Substituting the value of / to (A.l) gives 

|tr (LfcJd)l > di -h d 2 -h d3 - 2(d3 -h d 2 )(l - ql) > ds = ||Jd|| . 


□ 
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Appendix B. Proof of Claim 3.2 

Claim 3.2. Consider the equality RjRk+i = Fk- If we assume that the step length 
h is small enough such that the relative orientation RjRk+i between two adjacent 
time instances tk and tk+i is less than i.e., 

Il^fcll, IlCfcll < ^ where RjRk+i = e^^ = e^'= and Cfe,Cfc S 

then the aforementioned equality is satisfied if and only if the skew symmetric parts 
of both sides are identical. 


Proof. This claim can be proved using the Rodrigues’s formula. Let A := RjRk+i G 
SO(3) and B := Fk £ SO(3) then there exist vectors a, 6 £ R^ such that 

(B.l) ^^edlallSe) and B 

where Xe is a unit vector corresponding to the vector a; G R^. Using the Rodrigues’s 
formula [4] 

A = := / 3 x 3 + sin(||a||)de + dedj (cos(||a||) - 1) 

we obtain, the skew-symmetric parts of A and B as 

(B.2) S{A) := ^ ^ = sin(||a||)de and S{B) := ^ ^ = sin(||6||)&e. 

Given that S'(A) = S{B) and since we know that so (3) — R^, from (B.2) we have 
(B.3) sin(||a||)ae = sin(||&||)6e 


case 1: 

If ||6|| = 0 then from (B.3) sin(||a||)ae = 0 iff o = 0 which means a = b because 
||a|| < f ■ Using the exponential map (B.l) one can conclude that A = B. 

case 2: 


If ||6|| 7 ^ 0 then from (B.3) 


which is true iff 


sin(||6||)' 


= b. 


(B.4) 

(B.5) 


be = tte and 

sin(||Q||) ^ , 

sin(||6||) 


The map sin :] — §[—>■] — 1,1[ is bijective. Given that ||6|| < f, we conclude from 

(B.5) that sin(||a||) = sin(||6||). Hence 

(B.6) ||6|| = ||a|| 


From (B.4) and (B.6) it is clear that a = b. Using the exponential map (B.l), one 
can conclude that A = B. □ 


Appendix C. Proof of Theorem C.l 


Let us discuss a case in which, only momentum dynamics is considered . In this 
case, matching conditions (5.9) modifies to 

(C.l) C(y):=(Si Si ... Sfe Sfc ... Elv C^tm)^=0, 


where 

y := (nj 7 ([ •. • nj ... uj, 7 ^)^, s^+i := (n^+i - f^U k - huk )^, 
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^fc+l ■— ^ , Cnitm ■— ^ II _ ' 


We define the gradient of the matching conditions (C.l) as 


(C.2) 


where 


VyC{Y) := 


f-Ji 

0 

0 
0 

V 


/Cl 0 

-J2 IC 2 


Jk '■ = 




0 


Ik- 

hx3 


:= T^UoCmtm = 


0 0 

3x3 0 


0 

0 

— Jn-1 
0 
0 


0 \ 
0 


K-n-i 0 
—Jn /Cat 

0 Bf) 


JCk := 


^3x3 


"T 

k '~~‘k 


Pu 

:= Pn^Crntm = 


T) "T I 1 
‘^Ik^k . 


hx3 0 
0 0 


Theorem C.l. The Jacobian matrix VyC is non-singular if there exist a time 
instant k such that the control Uk is not saturated. 


Proof. The matrix "DyC is similar to the upper triangular matrix X>yC which is 
obtained using the elementary row transformations. 


f-Ji /Cl 0 
0 -J 2 IC 2 


(C.3) 


VyCiY) : = 


0 
0 
V 0 


0 

0 

—Jn-1 

0 

0 


0 \ 
0 


/Cjv -1 0 

—Jn /Cat 

0 X ) 


Note that the matrix VyCfY) is invertible if and only if Jk is invertible for k = 
1,2,..., N, and X is invertible. First we prove that Jk is invertible, 

(C.4) 


Jk = 


''Uk-i^k 

0 


'7k 

I: 


_ f ^Fk-i — hAfk-iFff_^Ilk-i^ 

3X3 J V 0 73x3 / 


Lemma C.2. The matrix {^Fk — hJ\fkF^J\.}^ is invertible if and only if the matrix 
(tr {FkJd) I 3 x 3 - FkJd) is invertible. 

A proof of Lemma C.2 is provided in Appendix D. By Lemma C.2 and Lemma 3.1, 
we know that the matrix (^Fk — hAfkF^Hk'^ is invertible. Hence, from (C.4) it is 
easy to conclude that the matrix Jk is invertible. Let us define the inverse of the 
matrix Jk as 

J-i := ^ fl3x3 

\ 0 I3X3 J \ Ck Uk 

where ^ 

ak — ^Fk hJJkFj^ H/,;^ , bk — hP-^^Uk, Ck — 

Taking the first order approximation in h, the matrix X can be approximated as 
X = Bf + B^J{^K.yJ 2^F2 ... Jff^lCN 


0 






ATTITUDE CONTROL WITH STATE AND CONTROL CONSTRAINTS 


23 


/3X3 0\ 

* sj ’ 

where s = such that = aNO-N-i ■ • ■ Ui- Now it is clear 

from (C.5) that the matrix X is invertible if the matrix s is invertible. As the matrix 
tti is invertible, so the matrix is invertible for all i = 0,1,. .. ,N. So, it is enough to 
prove that the matrix invertible for s being invertible. Note that 

the matrix bi is a negative semi-definite matrix for i = 1,2,... ^ N — 1. Assuming 
that there exist a time instant k such that the control Uk is not saturated means 
that bk is negative definite. So, it is concluded that the matrix E^o^ is 

negative definite and hence invertible. □ 

Appendix D. Proof of Lemma C.2 

Lemma C.2. The matrix (^Fk — hMkF^Yl}^ is invertible if and only if the matrix 
{tr (FkJd) hx 3 - FkJd) is invertible. 

Proof. We know that = F^hUkFk. Then using the implicit equation 

hHk = FkJd — JdF^, the matrix (^Fk — hMkF^ Ilfc^ can be rewritten as 

{Pk — hMkFj Ilfc^ = (^Fk — AfkFfl hUkFk^ 

= MkF,I (^Fk (A4)”^ - /Sfe) Fk 
= J^kF^ (tr {JdF^) /3x3 - JdF^ - wTfc) Fk 
= MkF^ (tr {FkJd) hxa — FkJd) Fk 

(D.l) = (tr {JdFj) /3x3 - JdFj)~" (tr (FkJd) hxs - FkJd) Fk. 

By (D.l) it is obvious to conclude that the matrix (tr {FkJd) hxa — FkJd) is invert¬ 
ible if and only if the matrix (^Fk — hMkF^U}^ is invertible. □ 
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